2  IFNg: Normalization, doublet discrimination, integration, clustering

2.1 Set up Seurat workspace

# Load libraries
library(data.table)
library(sctransform)
library(Seurat)
library(tidyverse)
library(miQC)   
library(SeuratWrappers)
library(fishpond)
library(Matrix)
library(patchwork)
library(alevinQC)
library(harmony)
library(scDblFinder)

# Set global options for Seurat v5 objects
options(Seurat.object.assay.version = 'v5')

2.2 Load previously saved merged object

merged.18279 <- readRDS("IFNg_scRNA_Part1.rds")

2.3 Normalize and scale data

Regress out mitochondrial contribution

merged.18279 <- NormalizeData(merged.18279)
merged.18279 <- FindVariableFeatures(merged.18279, 
                                    assay = "RNA", 
                                    layer = "data", 
                                    selection.method = "vst", 
                                    nfeatures = 5000)
merged.18279 <- ScaleData(merged.18279, vars.to.regress = "percent.mt")

2.4 Run PCA

merged.18279 <- RunPCA(merged.18279, npcs = 200)

PC_ 1 
Positive:  FABP5, PIM3, TNFRSF9, CD82, TNFRSF4, GNG4, PGAM1, PRDX1, ATP1B3, RAN 
       DUSP4, TNFRSF18, RPS17, NME1-NME2, C1QBP, MIR155HG, NDFIP1, NPM1, RPL7, RANBP1 
       ID2, SRM, HSP90AB1, HSPE1, PARK7, BACH2, PDCD1, HSPD1, REL, MAP2K3 
Negative:  IL32, CCL5, PLAAT4, LTB, TMSB4X, GZMA, CD74, PTPRC, PTPRCAP, CORO1A 
       IFITM1, MYL12A, CYBA, CYTIP, ANXA1, IFI6, GIMAP7, CD8A, RNF213, EVL 
       TMSB10, CD8B, LIMD2, SAMHD1, CTSW, ARHGDIB, GIMAP4, ISG15, STAT1, KLF6 
PC_ 2 
Positive:  RPL39, TPT1, RPS10, MIF, S100A10, LGALS1, S100A11, S100A6, TMSB10, PGK1 
       CD7, PLP2, TPI1, GPR183, PCSK1N, TMSB4X, LGALS3, ALDOA, LTB, LSR 
       RPL21P16, RPS10-NUDT3, CCR7, SPOCK2, CAPG, PLEKHA7, RPL36A, ENSG00000255508, UNC119, RPL21 
Negative:  CD38, ZEB2, ATXN1, CDK6, ELMO1, HLA-DRB5, CBLB, RABGAP1L, HLA-DRA, FUT8 
       LYST, PARP14, MIR4435-2HG, NEAT1, PRKCB, SIK3, PAM, RNF19A, RUNX2, ANKRD44 
       MBNL1, RBPJ, TIAM2, ENSG00000288632, ENTPD1, FYN, UTRN, SFMBT2, RNF213, TOX 
PC_ 3 
Positive:  MALAT1, BTG1, PNRC1, BNIP3L, REL, IL2, PPIAP22, TRIM22, PKIA, RPL21P75 
       RPL21P119, NINJ1, SLC2A3, LINC01619, RPL23AP42, PCED1B-AS1, ZBTB20, RPL21P16, MAML2, YPEL3 
       RPL7P9, RPL15P3, HNRNPA1P7, ENSG00000218426, BICDL1, RPL9P7, KLF2, IL7R, ENSG00000280441, P4HA1 
Negative:  ACTG1, TUBA1B, ACTB, TUBB, PSME2, CORO1A, MCM7, NME1-NME2, MCM5, GINS2 
       CALR, SLC25A5, DCTPP1, CCT6A, SRM, PA2G4, PRMT1, RAN, HSPE1, HSPD1 
       RANBP1, DYNLL1, GSTP1, SH3BGRL3, PCNA, DUT, LDHB, NHP2, MCM2, NPM1 
PC_ 4 
Positive:  SLC4A10, IFNG-AS1, PHLDA1, GNA15, IFNG, IL4I1, SLAMF1, HIF1A, SATB1, IL2RA 
       LGALS3, SLC2A3, RORA, PGK1, CD40LG, CYTIP, NFAT5, FURIN, BHLHE40, IL12RB2 
       DPP4, SLC16A3, ANKRD28, SYTL3, TNFAIP3, PFKFB3, NFKBIA, KLRB1, GPR183, IL2RB 
Negative:  IFITM1, RPS4Y1, ISG15, CTSW, CD8B, HLA-DPB1, LY6E, HLA-DPA1, HLA-DRB5, HLA-DRA 
       HLA-DRB1, PLAAT4, OASL, COTL1, CD27, KLRC4-KLRK1, SP110, IFIT3, IFIT1, C12orf57 
       HLA-DQA1, CD38, CCDC85B, OAS1, LIMD2, CD52, RPS10-NUDT3, IFIT2, HLA-DMA, ISG20 
PC_ 5 
Positive:  NOLC1, DDX21, NCL, HSPD1, HNRNPAB, EIF4A1, NOP16, FKBP4, EBNA1BP2, SRM 
       NPM1, TOMM40, C1QBP, SLC4A10, RPS17, PRMT1, NOP56, ODC1, ENSG00000249492, KCNQ5 
       CCT6A, SATB1, TRAP1, RPL17, DCTPP1, PA2G4, IFNG-AS1, CCT2, HSP90AB1, FEZ1 
Negative:  SH3BGRL3, C4orf3, S100A10, TPI1, ALDOA, BNIP3L, SIT1, LBH, GPI, IGFLR1 
       SRGN, CRIP1, GNA15, BNIP3, JUN, LGALS1, MIF, LSP1, CTSB, LDHA 
       LAG3, APOBEC3G, CD82, PGK1, NMB, TMSB4X, SLC16A3, BTG1, PGAM1, CD4 

ElbowPlot(merged.18279, ndims = 200, reduction = "pca")

2.5 Plot unintegrated UMAP

merged.18279 <- RunUMAP(merged.18279, 
                        reduction = "pca", 
                        dims = 1:50, 
                        reduction.name = "umap.unintegrated")

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session

15:24:48 UMAP embedding parameters a = 0.9922 b = 1.112

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by 'BiocGenerics'

15:24:48 Read 12057 rows and found 50 numeric columns

15:24:48 Using Annoy for neighbor search, n_neighbors = 30

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by 'BiocGenerics'

15:24:48 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

**************************************************|
15:24:52 Writing NN index file to temp file /tmp/RtmpzlkBZG/file37c5d823da5d58
15:24:52 Searching Annoy index using 1 thread, search_k = 3000
15:24:56 Annoy recall = 100%
15:24:57 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
15:25:00 Initializing from normalized Laplacian + noise (using RSpectra)
15:25:00 Commencing optimization for 200 epochs, with 532572 positive edges
15:25:09 Optimization finished

DimPlot(merged.18279, reduction = "umap.unintegrated", group.by = c("Site", "Patient", "Timepoint"))

2.6 Call preliminary clusters for the purposes of doublet discrimination

merged.18279 <- FindNeighbors(merged.18279, dims = 1:50, reduction = "pca")

Computing nearest neighbor graph

Computing SNN

merged.18279 <- FindClusters(merged.18279, 
                             resolution = 0.25, 
                             algorithm = 2)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 12057
Number of edges: 457888

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9125
Number of communities: 8
Elapsed time: 1 seconds

DimPlot(merged.18279, reduction = "umap.unintegrated", label = T)

2.7 Identify and remove doublets

This uses raw counts as input

2.7.1 Combine RNA layers

merged.18279[['RNA']] <- JoinLayers(merged.18279[['RNA']])

2.7.2 Convert seurat to sce and check colData

merged.18279.sce <- as.SingleCellExperiment(merged.18279, assay = "RNA")
merged.18279.sce

class: SingleCellExperiment 
dim: 61217 12057 
metadata(0):
assays(2): counts logcounts
rownames(61217): 5-8S-rRNA 5S-rRNA ... ZZEF1 ZZZ3
rowData names(0):
colnames(12057): P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG
  P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT ...
  P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC
  P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA
colData names(15): orig.ident nCount_RNA ... seurat_clusters ident
reducedDimNames(2): PCA UMAP.UNINTEGRATED
mainExpName: RNA
altExpNames(0):

colData(merged.18279.sce)

DataFrame with 12057 rows and 15 columns
                                                   orig.ident nCount_RNA
                                                  <character>  <numeric>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG   SeuratProject    20129.8
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT   SeuratProject    19535.9
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC   SeuratProject    19411.9
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT   SeuratProject    20429.8
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC   SeuratProject    19299.9
...                                                       ...        ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT SeuratProject       2048
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT SeuratProject       1499
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT SeuratProject       2483
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC SeuratProject       1177
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA SeuratProject       1471
                                                nFeature_RNA percent.mt
                                                   <integer>  <numeric>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG           4947    3.86697
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT           4384    2.45138
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC           4735    2.01548
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT           4958    2.64150
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC           4227    4.15930
...                                                      ...        ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT          965   4.736328
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT          778   1.667779
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT         1214   0.443012
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC          670   2.378929
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA          900   4.979606
                                                miQC.probability   miQC.keep
                                                       <numeric> <character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG          0.0756333        keep
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT          0.0434632        keep
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC          0.0573548        keep
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT          0.0424527        keep
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC          0.0975392        keep
...                                                          ...         ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT       0.05033876        keep
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT       0.02237223        keep
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT       0.27941627        keep
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC       0.00989129        keep
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA       0.07598866        keep
                                                    Patient        Site
                                                <character> <character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG          P108        IFNg
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT          P108        IFNg
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC          P108        IFNg
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT          P108        IFNg
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC          P108        IFNg
...                                                     ...         ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT        P104        IFNg
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT        P104        IFNg
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT        P104        IFNg
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC        P104        IFNg
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA        P104        IFNg
                                                  Timepoint   IpiCohort
                                                <character> <character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG       PostVax      5mgIpi
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT       PostVax      5mgIpi
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC       PostVax      5mgIpi
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT       PostVax      5mgIpi
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC       PostVax      5mgIpi
...                                                     ...         ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT     PostVax    2.5mgIpi
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT     PostVax    2.5mgIpi
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT     PostVax    2.5mgIpi
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC     PostVax    2.5mgIpi
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA     PostVax    2.5mgIpi
                                                      Assay
                                                <character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG           RNA
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT           RNA
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC           RNA
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT           RNA
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC           RNA
...                                                     ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT         RNA
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT         RNA
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT         RNA
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC         RNA
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA         RNA
                                                                Sample
                                                           <character>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG   P108_IFNg_PostVax_5m..
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT   P108_IFNg_PostVax_5m..
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC   P108_IFNg_PostVax_5m..
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT   P108_IFNg_PostVax_5m..
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC   P108_IFNg_PostVax_5m..
...                                                                ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT P104_IFNg_PostVax_2...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT P104_IFNg_PostVax_2...
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT P104_IFNg_PostVax_2...
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC P104_IFNg_PostVax_2...
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA P104_IFNg_PostVax_2...
                                                RNA_snn_res.0.25
                                                        <factor>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG                  0
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT                  0
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC                  0
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT                  0
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC                  0
...                                                          ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT                5
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT                1
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT                1
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC                1
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA                1
                                                seurat_clusters    ident
                                                       <factor> <factor>
P108_IFNg_PostVax_5mgIpi_RNA_ACTTACTTCCAAGCCG                 0        0
P108_IFNg_PostVax_5mgIpi_RNA_GCATACACAAAGGCGT                 0        0
P108_IFNg_PostVax_5mgIpi_RNA_TGCACCTAGACCCACC                 0        0
P108_IFNg_PostVax_5mgIpi_RNA_GGATGTTGTCTCTTAT                 0        0
P108_IFNg_PostVax_5mgIpi_RNA_TGCGGGTAGCATCATC                 0        0
...                                                         ...      ...
P104_IFNg_PostVax_2.5mgIpi_RNA_CCTAGCTGTCGTGGCT               5        5
P104_IFNg_PostVax_2.5mgIpi_RNA_CCCATACAGGTGATAT               1        1
P104_IFNg_PostVax_2.5mgIpi_RNA_ACTTTCATCACTTCAT               1        1
P104_IFNg_PostVax_2.5mgIpi_RNA_AGAGTGGGTTATTATC               1        1
P104_IFNg_PostVax_2.5mgIpi_RNA_TGGTTCCTCGTCTGAA               1        1

2.7.3 Run scDblFinder

Set the dbr.sd very high to better allow thresholds to be set based on misclassification rates per sample

merged.18279.sce <- scDblFinder(merged.18279.sce,
                                 samples = "Sample",
                 dbr.sd = 1,
                                 clusters = "seurat_clusters")

2.7.4 Inspect results

# Look at the classes
table(merged.18279.sce$seurat_clusters, merged.18279.sce$scDblFinder.class)

   
    singlet doublet
  0    2911     196
  1    2220     148
  2    2095     258
  3    1249     167
  4    1212      73
  5     751      11
  6     474      27
  7     253      12

table(merged.18279.sce$Sample, merged.18279.sce$scDblFinder.class)

                                
                                 singlet doublet
  P101_IFNg_PostVax_2.5mgIpi_RNA    2953     247
  P103_IFNg_PostVax_2.5mgIpi_RNA    1099     172
  P104_IFNg_PostVax_2.5mgIpi_RNA    3113     190
  P108_IFNg_PostVax_5mgIpi_RNA      4000     283

# Look at the scores
summary(merged.18279.sce$scDblFinder.score)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000525 0.0023453 0.0098138 0.1059070 0.0527471 0.9999471 

2.7.5 Save doublet classifications into main Seurat object

merged.18279$doublet_classification <- merged.18279.sce$scDblFinder.class

2.7.6 Count singlets and doublets

table(merged.18279$doublet_classification)

singlet doublet 
  11165     892 

table(merged.18279$doublet_classification, merged.18279$seurat_clusters)

         
             0    1    2    3    4    5    6    7
  singlet 2911 2220 2095 1249 1212  751  474  253
  doublet  196  148  258  167   73   11   27   12

2.7.7 Plot singlets/doublets in UMAP space

DimPlot(merged.18279,reduction = "umap.unintegrated", group.by = "doublet_classification")

2.7.8 Subset object to remove doublets and count remaining cells

merged.18279.singlets <- subset(merged.18279, doublet_classification %in% c("singlet"))
dim(merged.18279.singlets)

[1] 61217 11165

# Count remaining cells per initial cluster
table(merged.18279.singlets$seurat_clusters)

   0    1    2    3    4    5    6    7 
2911 2220 2095 1249 1212  751  474  253 

2.8 Remove cells with very high nCount_RNA values, set other final QC filters

merged.18279.singlets <- subset(merged.18279.singlets,
                                subset = nCount_RNA < 40000 & nCount_RNA > 3000 & nFeature_RNA > 1500)
dim(merged.18279.singlets)

[1] 61217 10320

2.9 Re-compute PCA

Re-scale data now that it has been subset

merged.18279.singlets[["RNA"]] <- split(merged.18279.singlets[["RNA"]], f = merged.18279.singlets$Sample)

merged.18279.singlets <- FindVariableFeatures(merged.18279.singlets, 
                                    assay = "RNA", 
                                    layer = "data", 
                                    selection.method = "vst", 
                                    nfeatures = 5000)
merged.18279.singlets <- ScaleData(merged.18279.singlets, vars.to.regress = "percent.mt")
merged.18279.singlets <- RunPCA(merged.18279.singlets, npcs = 200, verbose = TRUE)

2.10 Run Harmony integration

merged.18279.singlets <- IntegrateLayers(merged.18279.singlets, 
                                method = HarmonyIntegration, 
                                orig.reduction = "pca",
                                new.reduction = "integrated.harmony"
)

Warning: HarmonyMatrix is deprecated and will be removed in the future from the
API in the future

Warning: Warning: The parameters do_pca and npcs are deprecated. They will be ignored for this function call and please remove parameters do_pca and npcs and pass to harmony cell_embeddings directly.
This warning is displayed once per session.

Warning: Warning: The parameter tau is deprecated. It will be ignored for this function call and please remove parameter tau in future function calls. Advanced users can set value of parameter tau by using parameter .options and function harmony_options().
This warning is displayed once per session.

Warning: Warning: The parameter block.size is deprecated. It will be ignored for this function call and please remove parameter block.size in future function calls. Advanced users can set value of parameter block.size by using parameter .options and function harmony_options().
This warning is displayed once per session.

Warning: Warning: The parameter max.iter.harmony is replaced with parameter max_iter. It will be ignored for this function call and please use parameter max_iter in future function calls.
This warning is displayed once per session.

Warning: Warning: The parameter max.iter.cluster is deprecated. It will be ignored for this function call and please remove parameter max.iter.cluster in future function calls. Advanced users can set value of parameter max.iter.cluster by using parameter .options and function harmony_options().
This warning is displayed once per session.

Warning: Warning: The parameter epsilon.cluster is deprecated. It will be ignored for this function call and please remove parameter epsilon.cluster in future function calls. Advanced users can set value of parameter epsilon.cluster by using parameter .options and function harmony_options().
This warning is displayed once per session.

Warning: Warning: The parameter epsilon.harmony is deprecated. It will be ignored for this function call and please remove parameter epsilon.harmony in future function calls. If users want to control if harmony would stop early or not, use parameter early_stop. Advanced users can set value of parameter epsilon.harmony by using parameter .options and function harmony_options().
This warning is displayed once per session.

Transposing data matrix

Using automatic lambda estimation

Initializing state using k-means centroids initialization

Harmony 1/10

Harmony 2/10

Harmony converged after 2 iterations

2.11 Identify final clusters after integration

merged.18279.singlets <- FindNeighbors(merged.18279.singlets, dims = 1:30, reduction = "integrated.harmony")

Computing nearest neighbor graph

Computing SNN

merged.18279.singlets <- FindClusters(merged.18279.singlets, 
                             resolution = 0.7, 
                             algorithm = 2)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 10320
Number of edges: 352714

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8127
Number of communities: 11
Elapsed time: 1 seconds

merged.18279.singlets <- RunUMAP(merged.18279.singlets,
                        reduction = "integrated.harmony",
                        dims = 1:30,
                        reduction.name = "umap.harmony")

15:36:03 UMAP embedding parameters a = 0.9922 b = 1.112

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by 'BiocGenerics'

15:36:03 Read 10320 rows and found 30 numeric columns

15:36:03 Using Annoy for neighbor search, n_neighbors = 30

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by 'BiocGenerics'

15:36:03 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

**************************************************|
15:36:06 Writing NN index file to temp file /tmp/RtmpzlkBZG/file37c5d8437d41b
15:36:06 Searching Annoy index using 1 thread, search_k = 3000
15:36:10 Annoy recall = 100%
15:36:11 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
15:36:13 Initializing from normalized Laplacian + noise (using RSpectra)
15:36:13 Commencing optimization for 200 epochs, with 447378 positive edges
15:36:21 Optimization finished

table(merged.18279.singlets$seurat_clusters)

   0    1    2    3    4    5    6    7    8    9   10 
2023 1984 1845 1165 1093  778  470  373  355  130  104 

2.12 Plot clusters

DimPlot(merged.18279.singlets,reduction = "umap.harmony", label = TRUE)

DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "Patient")

DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "Site")

DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "Timepoint")

DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "IpiCohort")

DimPlot(merged.18279.singlets,reduction = "umap.harmony", group.by = "Sample") + NoLegend()

2.13 Save clustered object

saveRDS(merged.18279.singlets,"IFNg_scRNA_Part2.rds")

2.14 Get session info

sessionInfo()

R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Rocky Linux 8.10 (Green Obsidian)

Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so;  LAPACK version 3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/New_York
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] scDblFinder_1.14.0          SingleCellExperiment_1.22.0
 [3] SummarizedExperiment_1.30.2 Biobase_2.60.0             
 [5] GenomicRanges_1.52.1        GenomeInfoDb_1.36.4        
 [7] IRanges_2.34.1              S4Vectors_0.38.2           
 [9] BiocGenerics_0.46.0         MatrixGenerics_1.12.3      
[11] matrixStats_1.2.0           harmony_1.2.0              
[13] Rcpp_1.0.12                 alevinQC_1.16.1            
[15] patchwork_1.3.0             Matrix_1.6-4               
[17] fishpond_2.6.2              SeuratWrappers_0.3.19      
[19] miQC_1.8.0                  lubridate_1.9.3            
[21] forcats_1.0.0               stringr_1.5.1              
[23] dplyr_1.1.4                 purrr_1.0.2                
[25] readr_2.1.5                 tidyr_1.3.1                
[27] tibble_3.2.1                ggplot2_3.4.4              
[29] tidyverse_2.0.0             Seurat_5.1.0               
[31] SeuratObject_5.0.2          sp_2.1-3                   
[33] sctransform_0.4.1           data.table_1.15.0          

loaded via a namespace (and not attached):
  [1] spatstat.sparse_3.0-3     bitops_1.0-7             
  [3] httr_1.4.7                RColorBrewer_1.1-3       
  [5] tools_4.3.1               utf8_1.2.4               
  [7] R6_2.5.1                  DT_0.31                  
  [9] lazyeval_0.2.2            uwot_0.1.16              
 [11] withr_3.0.0               GGally_2.2.1             
 [13] gridExtra_2.3             progressr_0.14.0         
 [15] cli_3.6.2                 spatstat.explore_3.2-6   
 [17] fastDummies_1.7.3         labeling_0.4.3           
 [19] spatstat.data_3.0-4       ggridges_0.5.6           
 [21] pbapply_1.7-2             Rsamtools_2.16.0         
 [23] R.utils_2.12.3            scater_1.28.0            
 [25] parallelly_1.37.0         limma_3.56.2             
 [27] generics_0.1.3            BiocIO_1.10.0            
 [29] gtools_3.9.5              ica_1.0-3                
 [31] spatstat.random_3.2-2     ggbeeswarm_0.7.2         
 [33] fansi_1.0.6               abind_1.4-5              
 [35] R.methodsS3_1.8.2         lifecycle_1.0.4          
 [37] yaml_2.3.8                edgeR_3.42.4             
 [39] Rtsne_0.17                grid_4.3.1               
 [41] promises_1.2.1            dqrng_0.3.2              
 [43] crayon_1.5.2              shinydashboard_0.7.2     
 [45] miniUI_0.1.1.1            lattice_0.22-5           
 [47] beachmat_2.16.0           cowplot_1.1.3            
 [49] pillar_1.9.0              knitr_1.45               
 [51] metapod_1.8.0             rjson_0.2.21             
 [53] xgboost_1.7.7.1           future.apply_1.11.1      
 [55] codetools_0.2-19          leiden_0.4.3.1           
 [57] glue_1.7.0                remotes_2.4.2.1          
 [59] vctrs_0.6.5               png_0.1-8                
 [61] spam_2.10-0               gtable_0.3.4             
 [63] xfun_0.42                 S4Arrays_1.2.0           
 [65] mime_0.12                 survival_3.5-8           
 [67] statmod_1.5.0             bluster_1.10.0           
 [69] ellipsis_0.3.2            fitdistrplus_1.1-11      
 [71] ROCR_1.0-11               nlme_3.1-164             
 [73] RcppAnnoy_0.0.22          irlba_2.3.5.1            
 [75] vipor_0.4.7               KernSmooth_2.23-22       
 [77] colorspace_2.1-0          nnet_7.3-19              
 [79] tidyselect_1.2.0          compiler_4.3.1           
 [81] BiocNeighbors_1.18.0      DelayedArray_0.26.7      
 [83] plotly_4.10.4             rtracklayer_1.60.1       
 [85] scales_1.3.0              lmtest_0.9-40            
 [87] digest_0.6.34             goftest_1.2-3            
 [89] spatstat.utils_3.0-4      rmarkdown_2.25           
 [91] RhpcBLASctl_0.23-42       XVector_0.40.0           
 [93] htmltools_0.5.7           pkgconfig_2.0.3          
 [95] sparseMatrixStats_1.12.2  fastmap_1.1.1            
 [97] rlang_1.1.3               htmlwidgets_1.6.4        
 [99] shiny_1.8.0               DelayedMatrixStats_1.22.6
[101] farver_2.1.1              zoo_1.8-12               
[103] jsonlite_1.8.8            BiocParallel_1.34.2      
[105] R.oo_1.26.0               BiocSingular_1.16.0      
[107] RCurl_1.98-1.14           magrittr_2.0.3           
[109] modeltools_0.2-23         scuttle_1.10.3           
[111] GenomeInfoDbData_1.2.10   dotCall64_1.1-1          
[113] munsell_0.5.0             viridis_0.6.5            
[115] reticulate_1.35.0         stringi_1.8.3            
[117] zlibbioc_1.46.0           MASS_7.3-60.0.1          
[119] plyr_1.8.9                flexmix_2.3-19           
[121] ggstats_0.5.1             parallel_4.3.1           
[123] listenv_0.9.1             ggrepel_0.9.5            
[125] deldir_2.0-2              Biostrings_2.68.1        
[127] splines_4.3.1             tensor_1.5               
[129] hms_1.1.3                 locfit_1.5-9.8           
[131] igraph_2.0.2              spatstat.geom_3.2-8      
[133] RcppHNSW_0.6.0            reshape2_1.4.4           
[135] ScaledMatrix_1.8.1        XML_3.99-0.16.1          
[137] evaluate_0.23             scran_1.28.2             
[139] BiocManager_1.30.22       tzdb_0.4.0               
[141] httpuv_1.6.14             RANN_2.6.1               
[143] polyclip_1.10-6           future_1.33.1            
[145] scattermore_1.2           rsvd_1.0.5               
[147] xtable_1.8-4              restfulr_0.0.15          
[149] svMisc_1.2.3              RSpectra_0.16-1          
[151] later_1.3.2               viridisLite_0.4.2        
[153] beeswarm_0.4.0            GenomicAlignments_1.36.0 
[155] tximport_1.28.0           cluster_2.1.6            
[157] timechange_0.3.0          globals_0.16.2